Structural phase transition and orientation-strain glass formation in anisotropic 
particle systems with impurities in two dimensions 
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Using a modified Lennard- Jones model for elliptic particles and spherical impurities, we present 
results of molecular dynamics simulation in two dimensions. In one-component systems of elliptic 
particles, we find an orientation phase transition on a hexagonal lattice as the temperature T is low- 
ered. It is also a structural one because of spontaneous strain. At low T, there arise three martensitic 
variants due to the underlying lattice, leading to a shape memory effect without dislocation forma- 
tion. Thermal hysteresis, a minimum of the shear modulus, and a maximum of the specific heat 
are also found with varying T. With increasing the composition c of impurities, the three kinds of 
orientation domains are finely divided, yielding orientation-strain glass with mesoscopically ordered 
regions still surviving. If the impurities are large and repulsive, planar anchoring of the elliptic par- 
ticles occurs around the impurity surfaces. If they are small and attractive, homeotropic anchoring 
occurs. Clustering of impurities is conspicuous. With increasing the anchoring power and/or the 
composition of the impurities, positional disorder can also be enhanced. We also investigate the 
rotational dynamics of the molecular orientations. 



I. INTRODUCTION 

Certain anisotropic molecules such as N2, Ceo, and 
KCN form a cubic crystal and, at lower tempera- 
tures, they undergo an orientation phase transition 
with a specific-heat peakpEl, where the crystal struc- 
ture changes to a noncubic one. Furthermore, mix- 
tures of anisotropic particles^ such as (KCN) a: (KBr)i_ a: 
and one-component systems of globular molecules^ such 
as ethanol and cyclohcxanol become orientation glass. 
In such glass, the phase ordering should occur only 
on small spatial scales with mesoscopically heteroge- 
neous orientation fluctuations. Because of anisotropic 
molecular shapes, there should be a direct (proper) cou- 
pling between the molecular orientations and the lat- 
tice deformations®. In fact, the shear modulus becomes 
small ar ound the orientational order-disorder or glass 
transitiorpEl. These systems thus exhibit singular acous- 
tic and plastic behaviors^, but there has been no system- 
atic experiment in the nonlinear response regime. Many 
of these anisotropic molecules have dipolar moments also, 
yielding dielectric anomaly near the transition. 

In metallic alloys, a structural phase transition arises 
from the displacements of the atoms in each unit cell 
from their equilibrium positions in the high-symmetry 
phase. Some alloys undergo a martensitic phase transi- 
tion gradually from a high-temperature phase to a low- 
temperature phase over a rather wide temperature range 
and, at sufficiently low temperatures, they are com- 
posed of multiple martensitic variants or domains-^ESD. 
In particular, a system of off-stoichiometric intermetallic 
Tiso-zNiso+a; has been studied extensively^M Even at 
x = 1.5, it becomes strain glass, exhibiting the shape- 
memory effect and the superclasticity, where strain het- 
erogeneities with sizes of order 10 nm were observed^. 
As a similar example, metallic ferroelectric glass, called 
relaxor, exhibits large dielectric response to applied elec- 
tric fielcP^HS] where the electric polarization and the 



lattice deformations are coupled and frozen polar nan- 
odomains are produced in the presence of the composi- 
tional disorder in the perovskite structure. 

In soft matter, impurities often strongly disturb or in- 
fluence phase transitions. Examples of impurities are 
filler particles in phase-separating polymer blends^, mi- 
croemulsions in nematic liquid crystals^, and crosslink 
irregularities in polymer gela^HSSI. p or ggi S; some au- 
thors developed random crosslink model*^. Moreover, 
in gels with liquid crystal solvent ^211221, ^he isotropic- 
nematic phase transition is analogous to the orientation 
phase transition in solids, where the coupling between 
the molecular orientation and the elasticity leads to sin- 
gular elastic behavior. In such liquid crystal gels, ne- 
matic polydomains are produced by random crosslinkage 
and polydomain-monodomain transitions are induced by 
applied stress or electric fielcPH as numerically studied 
by Uchida 24 using quenched random stress. The poly- 
domains obviously correspond to the mesoscopic orienta- 
tion heterogeneities in solids. We also mention exper- 
iments of crystal formation and glass transition using 
elongated colloidal particles in three dimensions^ and 
in two dimension^. 



The mesoscopic heterogeneities produced by impu- 
rities are widely recognized in various solid and soft 
materials. We may mention two previous approaches. 
One is based on a ra ndom field coupled to the or- 
der parameter^E^IMIZHSD j n particular, Vasseur and 
Lookmarl22i introduced a spin glass theory supplemented 
with the elastic interaction (the long-range interaction 
among the order parameter ip mediated by the elastic 
deformations}^. The other is a phase-field (Ginzburg- 
Landau) theory with a random critical temperature and 
the elastic interaction^"^!, where the quadratic term 
(oc ip 2 ) in the free energy has a random coefficient. In 
these theories, the impurities are governed by an arti- 
ficial random distribution without spatial correlations. 
Therefore, they lack microscopic physical pictures of the 



impurity disordering. From our viewpoint, microscopic 
approaches are particularly needed when each impurity 
strongly perturbs the local order parameter. 

To perform first-principle calculations of the meso- 
scopic heterogeneities, we start with Lennard- Jones sys- 
tems composed of anisotropic host particles and impuri- 
ties to create orientationally disordered and ordered crys- 
tal states, (i) In such states we may examine the degree 
of heterogeneities by changing the impurity composition, 
(ii) We may describe a tendency of impurity clustering 
or aggregation 33 , which depends on the cooling rate from 
liquid. It apparently governs the degree of vitrification, 
for example, in water containing a considerable amount of 
salipH (iii) We also note that the positional disorder and 
the orientation disorder have been discussed separately in 
the literature. In this paper, they appear simultaneously, 
though the former is weaker than the latter. 

The organization of this paper is as follows. In Sec. II, 
we will present the backgrounds of our theory and sim- 
ulation. In Sec. Ill, we will give simulation results for 
one-component systems of elliptic particles forming crys- 
tal to examine the orientation transition. In Sec. IV, we 
will treat mixtures of elliptic particles and larger repul- 
sive impurities, where the elliptic particles are aligned in 
the planar alignment around the impurities^. In Sec.V, 
we will examine the orientation dynamics of the elliptic 
particles. In Sec. VI, we will treat small attractive impu- 
rities, which tend to form aggregates and solvate several 
elliptic particles in the homeotropic alignment!^. 



II. THEORETICAL AND SIMULATION 
BACKGROUNDS 

We propose a simple microscopic model of binary mix- 
tures in two dimensions, which exhibits orientation phase 
transitions and glass behavior. We do not introduce the 
dipolar interaction supposing nonpolar molecules. 



A. Angle-dependent potential 

In our model, the first and second components are 
composed of elliptic and spherical particles, respectively. 
Their numbers are iVi and N2, where N = Ni + N 2 — 
4096 in this paper. The composition is defined by 



N 2 /N, 



(1) 



which is either of 0, 0.05, 0.1, 0.15, 0.2, or 0.3 in this 
paper. Thus the particles of the second component con- 
stitute impurities. The particle positions are written as 
Ti (i = 1, • • • , N). The orientation vectors of the elliptic 
particles may be expressed in terms of angles 6i as 



where i = 1, 



n, : = (cos^sin^), 
,JVi. 



(2) 



The pair potential [Zy between particles i £ a and 
j € j3 (a,j3 = 1,2) is a truncated modified Lennard- 
Jones potential. That is, for r,-j > r c = ioi it is zero, 



while for ry 


< r c = 3<7i it reads 




Uij = 4e 


(i + ^.)5£-(i + s<i)5r 


-Cir (3) 


L - y • ij j 


Here, r$ — r 


j = TijTij with Vij = |ry|. In terms of the 



diameters <j\ and 02 of the two species, we define 

o a /3 = (o a +o p )/2. (4) 

In Eq(3), Cy is the value of the first term at r = r c , 
ensuring the continuity of £/y . The e is the characteristic 
interaction energy. 

The particle anisotropy is taken into account by the 
angle factors Ay and -By, which depend on the relative 
direction f y = r~- l rij and the orientations n^ and rij of 
the elliptic particles. There can be a variety of their forms 
depending on the nature of the anisotropic interactions. 
Throughout this paper, we assume the following form, 



Aij = x5ai{rii ■ Tijf + x8pi(nj ■ rij) 2 



(5) 



where \ is the anisotropy strength of repulsion. The 8 a \ 
((5/31) is equal to 1 for a = 1 (/3 = 1) and for a = 2 
(,3 = 2). Thus, in the right hand side, the first (second) 
term is nonvanishing only when i (j) belongs to the first 
species. In Sec. IV, we treat large spherical impurities 
repelling the elliptic particles by setting <72/<7i > 1 and 
Bij =0. If x > in this case, there appears a tendency of 
parallel alignment of the elliptic particles at the impurity 
surfaces 35 . On the other hand, in Sec. VI, we assume 
oil a \ < 1 and 

B^ = (5 a iS/32(ni ■ rij) 2 + C^a2^i("j • r^) 2 , (6) 

where £ is the anisotropy strength of attraction. In this 
case, the attractive interaction is anisotropic only be- 
tween the elliptic particles and small spherical impurities 
and, if £ > 0, there appears a tendency of homeotropic 
alignmenlpS at the impurity surfaces. 

The total energy is written as T-L = K + U, where U is 
the potential energy and K is the kinetic energy, 



l<i<j<N 



h 



%\ 2 , 



(7) 
(8) 



Ki<N 



l<i<N! 



where f^ = dri/dt, 9i = d$i/dt, mi and m,2 are the 
masses, and I\ is the moment of inertia of the first com- 
ponent. In this paper, we set mi = m 2 = m. The 
Newton equations of motion are now written as 



_d_ 
dr, 
d 



K 



_d_ 
dr, 



U, 



n"i 



09,, 



K 



00, 



u, 



(9) 

(10) 



where i £ a, r t = d 2 ri/dt 2 , and Qi = d 2 9i/dt 2 . The sec- 
ond line holds for the first component (i — 1, ■ • • , Ni). 
However, since we treat equilibrium or at least nearly 
steady states, we attach a Nose-Hoover thermostat)^ 
to all the particles by adding the thermostat terms in 
Eqs.(9) and (10). Unless confusion may occur, space, 
time, and temperature will be measured in units of o\ , 



To = 0"! 



y/mx/e, 



(11) 



and e/kg, respectively, where ks is the Boltzmann con- 
stant. Stress (and elastic moduli) will be measured in 
units of e/a 2 . 

From Eqs.(3), (5), and (6) the elliptic particles have 
angle-dependent diameters. Let the particles i and j be- 
long to the first species. Then minimization of E/y in 
Eq.(3) with respect to r tj gives ry = 2 1 / 6 (1 + Ay) 1 / 6 ^. 
Thus the shortest diameter a s is given for the perpen- 
dicular orientations (n$ • ry = n» ■ ry = 0), while the 
longest diameter at by the parallel orientations (rij • ry = 



±1) so that 



= 2 1 / 6 



o\. 



(if 



= (l + 2 X ) 1/6 2 1/( 



v\- 



(12) 



The ratio of these lengths (the aspect ratio) is given by 
ai/a s = (1 + 2x) 1 ^ 6 - For example, ai/a s is equal to 1.14 
for x = 0-6 and to 1.23 for x = 1.2. We estimate the 
effective molecular area and the momentum of inertia of 
the elliptic particles as 

Si = TTa s a £ /4, h = (aj + a^rra/i. (13) 

In this paper, we fix the overall packing fraction as 

P ack = (JViSi + N 2 S 2 )/V = 0.95, (14) 

where S 2 — 7r2 1//3 (r|/4 and V is the system volume. Then 
the system length is L — V x l 2 = 70cti. 

Our potential is analogous to the Gay-Berne potential 
for anisotropic molecules used to simulate mesophases 
of liquid crystals^ and the Shintani-Tanaka potential 
with five-fold symmetry used to study frustrated parti- 
cle configurations at high densities®. It is worth noting 
that angle-depende nt pot entials have been used for lipids 
forming membranes ! 39 ! 40 ! 



B. Coarse-grained orientation order parameter 

For each particle % of the first species (i = 1, • • • , N\), 

we may introduce the orientation tensor Q i = {Qi^ u } 
(//, v = x,y) in terms of the orientation vectors n^ as 

Q i = (1 + n^in.n, + ^ rijUj)- / /2 

j £ bonded 



qi(didi- I /2), 



(15) 



where I = {Suu} is the unit tensor and dj is the direc- 
tor with |cL I = 1. The summation is over the bonded 



particles (|ry| < 3cti) of the first species with n l h being 
the number of these particles of order 20. If a hexagonal 
lattice is formed, it includes the second nearest neighbor 
particles. The angle of the director is defined by 

d t = (cos ipi, sin. <pi), (16) 

in the range < ifi < w. The amplitude <& is given by 

q? = 2j2Qlv (17) 



\.l,V 



We will calculate the average over the elliptic particles, 



l<J<iV! 



(18) 



which represents the overall degree of orientation order. 
The angle ifi varies more smoothly than 9i, but they 
coincide in ordered domains at low T. The q 2 is of order 
0.1 in disordered states due to the thermal fluctuations, 
but it increases up to unity within domains at low T. 
As a merit in visualization, q 2 is small in the interface 

regions at low T (see the right panels of Fig.l). 
<-> 
Since the tensor Q i is symmetric and traceless, its 
components are written as Qi XX = —Qi yy = Q«2/2 and 
Qixy = Qiyx = Qi3/2. In terms of ifi we have 



Q12 = qt cos(2(pi), Q a = q t sin(2ipi). 



(19) 



These variables change with respect to a rotation of the 
reference frame by an angle tp as^Sl 



Q'i2 = 
Qi3 = 



) l2 cos(2^) 
) i3 cos{2tp) 



} l3 sm(2ip), 
)i2sin(2i/0- 



We also introduce the following density variables as 



iei 



(20) 

(21) 
(22) 



iei 



We will calculate the following structure factor, 

S Q (k) = (\Q 2k \ 2 ) 

= (IQ 3 fel 2 )' 



(23) 



where Q 2 fc and Q^ are the Fourier components of Q 2 (r) 
and Qs(r), respectively. From Eq.(20) the structure 
factor of Q 2 and that of Q3 coincide under the ro- 
tational invariance of the system (without stretching), 
leading to the second line of Eq.(23). If there is no 
anisotropic overall strain, the isotropy holds for k much 
larger than the inverse system length, leading to Eq.(23) 
and (Q 2k Q* 3k ) = 0. 



8=Ti<^ C=0,X=0.6 



T=0.09 




FIG. 1. Orientation angles di in the range < 8i < it (left) 
and order parameter amplitudes qf (right) of all the parti- 
cles on a lattice in the xy plane with c = and x — 0.6 
for T = 0.09, 0.074, 0.07, and 0.04 from above. As T is 
lowered, orientation order develops gradually with lattice de- 
formations. 



III. ORIENTATION PHASE TRANSITION IN 
ONE-COMPONENT SYSTEMS 
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FIG. 2. Left: Expanded snapshot of 0i around a junc- 
tion point of three variants in the bottom panel of Fig.l at 
T = 0.04. The angles among the three lines are 65, 145, and 
150 degrees, being approximately multiples of n/6. Right: 
Hexagonal lattice structure in an ordered variant composed 
of isosceles triangles for \ = 0.6 and T = 0.04. 
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FIG. 3. Left: Structure factor of the orientation fluctua- 
tions Sg(k) in Eq. (23) vs k with c = and \ = 0.6 for 
T = 0.07, 0.072, 0.074, 0.076, 0.08, and 0.09, which grows alge- 
braically for small k < 0.3 in the BKT phase. Right: Average 
amplitude (g 2 ) in Eq.(18) vs T for x = 0.6, 1.2, 1.8, and 2.4, 
which increases gradually but steeply at low T. 



disordered and ordered crystal phases. Solids in the ori- 
entati onally d isordered phase have been called "plastic 
solids'^^EIHMl To understand singular mechanical prop- 
erties of TiNi around its martensitic phase transition, 
Ding et al. performed molecular dynamics simulation 
on mixtures of two species of spherical particles. 



In this section, we treat pure (one-component) sys- 
tems of the elliptic particles (c = 0). We assume not 
large values of x(< 2.4) such that the crystallization first 
occurs at T = T m ~ 1 with random molecular orien- 
tations. Far below T m , we study an orientation phase 
transition on a hexagonal lattice and singular mechani- 
cal behavior specific to multi- variant states. A number of 
author a 41 " 43 ! numerically examined the phase behavior of 
one-component hard rod systems in three dimensions in 
the plane of the aspect ratio and the density. If the parti- 
cles are rather close to spheres, they found orientationally 



A. Variant formation and 
Berezinskii-Kosterlitz-Thouless phase 

In Figs. 1-3, we show our simulation results at fixed 
volume under the periodic boundary condition. Assum- 
ing a Nose-Hoover thermostalpSl, we started with a liquid 
at T = 2, quenched the system to T = 0.35 below the 
melting, and annealed it for 9000to. We then lowered 
T to a final low temperature. Here, even if the cooling 
rate was varied after the crystal formation (in the range 
T < 0.35), essentially the same results followed. That is, 



there was no history-dependent behavior. 

In Fig.l, we show the orientation angles 0, of all the 
particles in the range < Oi < ir (left) and the or- 
der parameter amplitudes qf in Eq.(17) (right) at T — 
0.09, 0.074, 0.07, and 0.04. With lowering T, three equiv- 
alent variants emerge due to the underlying hexagonal 
lattice. Their areal fractions are all nearly equal to 1/3. 
For T = 0.074 the time scale of the patterns is of order 
10 4 , while for T = 0.07 and 0.04 the patterns are frozen 
even on time scale of 10 5 . At low T, the junction angles, 
at which two or more domain boundaries intersect, are 
multiples of n/6. As illustrated in Fig. 2, this geometrical 
constraint arises from the orientation-lattice coupling. It 
serves to pin the domain growth at a characteristic size 
even without impurities 20 . Similar pinned domain pat- 
terns have been observed on hexagonal planes^ and were 
reproduced by phase-field simulation 4 ^ We may define 
the surface tension 7 on the interfaces far from the junc- 
tion regions. In our model, 7 ~ 0.1e/er 2 for \ = 0.6 and 
7 ~ 0.2e/cr 2 for * = 1.2 at low T. 

In Fig. 3, the structure factor Sq(k) in Eq.(23) vs k 
and the average (q 2 ) in Eq.(18) vs T are displayed for 
X = 0.6,1.2,1.8, and 2.4. Here, the orientation order 
develops continuously in a narrow temperature range, 






T 2 ( X )<T<T 1 (x), 



(24) 



where T 2 ~ 0.070 and T x ~ 0.076 for X = 0.6. In our sim- 
ulation, T\ and T 2 increase with increasing \- They are 
determined as crossover temperatures. In this temper- 
ature window, a Berezinskii-Kosterlitz-Thouless (BKT) 
phas e 47 ! 48 ! is realized between the low-temperature or- 
dered phase for T < T 2 and the high-temperature disor- 
dered phase for T > T\, where the orientation fluctua- 
tions are strongly enhanced at long wavelengths. In our 
model, each elliptic particle on a lattice point behaves 
as a rotator in the XY spin model under a symmetry- 
breaking free energy AF — — ^\ h p cos(p0i) with p = 6, 
which arises from the under lying crystal structure 4 ^. In 
accord with the theory^lMSJ the structure factor Sq(k) 
in Eq.(23) grows algebraically as 



S Q (k) ~ k- 



(k < 0.5) 



(25) 



in the temperature range (24) with 77 depending on T 
(r) = 0.05 at T — 0.074). As regards dynamics, the orien- 
tation fluctuations migrate in space on rather rapid time 
scales slightly below T\, but are frozen for T < T 2 (see 
Fig. 13 below). Considerably below T 2 , the three variants 
become distinct with sharp interfaces. Previously, for 
two-dimensional hard rods, Bates and Frenkel 49 found 
a Kosterlitz-Thouless phase transition between the ne- 
matic phase and the isotropic phase for large aspect ra- 
tios and for low densities. 

As illustrated in the right panel of Fig. 2, the orienta- 
tion order induces lattice deformations. In ordered states 
at low T, each variant is composed of isosceles triangles 
elongated along its orientated direction parallel to one of 
the crystal axes. At low T, their side lengths b and c 




0.09 



C.09 



FIG. 4. Order parameter amplitude (q 2 ) (top) and strain 
e (bottom) under fixed applied stress a a = +0, 0.05, 0.10, 
and 0.15 in units of s/a 2 , where c = and x — 0.6. The 
temperature T was first decreased from 0.1 to 0.02 and it 
was then increased back to 0.1, where dT/dt = =p4 x 10 -6 . 
Hysteretic behavior appears between the cooling and heating 
paths. 



are (6, c) = (1.21,1.11) for X = 0.6 and (1.27,1,11) for 
X = 1.2 under the periodic boundary condition, while we 
have (b,c) = (1.28, 1.12) at zero stress. Thus this orien- 
tation transition is also a structural or martensitic one 
with spontaneous lattice deformations. 



B. Mechanical properties and specific heat for c = 

We have also performed simulation at a fixed stress*"", 
which allows an anisotropic shape change of the sys- 
tem at a structural phase transition. In Figs. 4-7, we 
assumed a Nose-Hoover thermostatP3 and a Parrinello- 
Rahman barostat^ under the periodic boundary condi- 
tion. Namely, we controlled the temperature T and the 
stress along the y axis written as 



(26) 



where (• • •) represents the space average. The y axis is 
taken to be in the perpendicular direction in the figures. 
Hereafter, a a will be measured in units of e/af. When 
a a was held fixed at a positive value for a long time, a 
single- variant state elongated along the y axis was even- 
tually realized at low T. This was the case even for very 




FIG. 5. Young's modulus E e in Eq.(28) divided by 4 (left) 
and the isobaric specific heat C p in Eq.(32) (right) for c = 
on the cooling and heating paths in Fig. 4 in the nearly stress- 
free condition (a a = 10 -3 ). Here E e /4 is nearly equal to the 
effective shear modulus /x e in Eq.(29). Softening against shear 
deformations and large energy fluctuations are conspicuous at 
the orientation transition. 
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FIG. 6. Shape memory effect under uniaxial deformations 
along the y axis without impurities (c = 0), where T = 0.02 
and x — 0.6. Left: Strain e vs applied stress a a in units of 
e/crf. For a a > 0.075, there remains only the variant elon- 
gated along the y axis. After this cycle, the residual strain 
vanishes upon heating to T = O.f. Right: Fractions of the 
three variants during the cycle, which are stretched along the 
three crystal axes. 



small positive a a (~ 10~ 3 ), since it serves as a symmetry- 
breaking field. We also carried out many simulation runs 
exactly setting a a = 0, where a few domains often re- 
mained in the final state depending on the initial condi- 
tions (not shown in this paper). When a a is controlled, 
the system length L y along the y axis should be calcu- 
lated. The strain e is defined as 



Ly/L 



yo 



1, 



(27) 



where L y o is a reference system length to be specified 
below. We may define Young's modulus by 



E P 



l/(de/da a ) 



(28) 



even in the nonlinear regime. Note that Young's modulus 
is written as E — 4Kfi/(K + fi) in the linear regime in 
terms of the bulk modulus K and the shear modulus fi 



in two dimensions. We may introduce the effective shear 
modulus \i e replacing E by E e as 



^ = E e /{4 - E e /K) £* E e /4, 



(29) 



where we have assumed E e <C K. 

Substantial thermal hysteresis during cooling and heat- 
ing has been observed in alloys around martensitic phase 
transitions^]. Ding et al. also found thermal hysteresis 
numerically . In Fig. 4, we show thermal hysteresis in 
our system for c = 0. That is, fixing a a , we decreased 
T from 0.1 to 0.02 with a very slow cooling rate given 
by dT/dt = —4 x 10~ 6 , where the variant elongated 
along the y axis became dominant at low T. We then 
increased T back to the initial high temperature with 
dT/dt = 4 x 10~ 6 . The curves of a a = +0 arc those 
with a small symmetry breaking stress (= 10~ 3 ). The 
reference length L y0 in Eq.(27) is that at T = 0.09 equal 
to 73.0, 73.5, 75.7, and 76.2 for a a = +0, 0.05, 0.1, and 
0.15, respectively. Hysteretic behavior can be seen in the 
degree of orientation (q 2 ) and the strain e. The width of 
the hysteresis loop is maximum for a a = +0 and shrinks 
to vanish for a a > 0.15. The transition at a a — +0 be- 
tween the orientationally disordered and ordered states 
is shifted to lower temperatures by 0.01 than in the fixed- 
volume case in Fig.l. See Remark (4) in Sec. VII for dis- 
cussions on the stability of quasi-equilibrium states in 
Fig.4. 

When (q 2 ) is small, we may use the linear elasticity 
relations in two dimensions, 

K[ei - a(T - T Q )] + ^{2e - ei) = a a , 

K[ei -a(T- T )] -n(2e-ei) = 0, (30) 

where e\ is the dilation strain, a is the thermal expansion 
coefficient, and Tq is a reference temperature . The small 
slope of the curves of e at a a = +0 in the disordered 
regions in Fig.4 arise from the thermal expansion. From 
these relations we obtain 



e l = a a /2K + a(T-T ), 

e = (K + /iK/4A> + a(T - T )/2 



(31) 



The data in Fig.4 yield K = 20, a = 0.6, and \x = 2 with 
T = 0.09 in the disordered phase. 

In the left panel of Fig. 5, we show Young's modulus E e 
in Eq.(28) on the cooling and heating paths of a a = +0 
in Fig.4. To calculate it, we superimposed a small stress 
(= 10~ 2 ) to the much smaller symmetry-breaking stress 
(= 10~ 3 ). Remark ably, E e becomes very small around 
T = 0.06 for c = d 2 ° l 51 l Similar minimum behavior of 
the shear modulus has been observed near the orientation 
and glass transition o 1 ! ' 5 \ where the minimum depends on 
the mixture composition. Previously, using the correla- 
tion function expression, Murat and Kantor^S calculated 
the elastic constant to find its softening toward the ori- 
entation transition in two-dimensional ellipsoid systems. 
Nonlinear response behavior should appear even for very 
small applied strains near the transition. Additionally, in 
the right panel of Fig. 5, we display the isobaric specific 



heat in the nearly stress-free condition (a a — 10 3 ) along 
the cooling and heating paths expressed as 



9 ^H" X=1 2, a 2 / 0l - 1.2, T=0.05 

c=0.05 



C P = ((H-(H)) 2 )/VT 2 



(32) 



where H = K + U is the total energy (see Eqs.(7) and 
(8)) and SW = H — (H) is its deviation. It is peaked at 
T = 0.06 indicating enhancement of the energy fluctua- 
tions at the transition. Such specific heat anomaly has 
been measured near the orientation transition^. We also 
calculated the constant- volume specific heat Cy using the 
data in Fig.l to find a similar peak around T = 0.073 (not 
shown in this paper). 

Next, we illustrate the shape memory effect taking 
place without dislocations. In Fig. 6, we increased a a 
from to 0.1 and then decreased a a back to at T = 0.02, 
where da a /dt = ±7 x 10~ 6 with + being on the stretch- 
ing path and — being on the return path. In this slow 
cycle, the system remained in quasi-static states. In the 
definition of e in Eq.(27), L y0 is the initial system length 
(= 72). At t = 0, the fractions of the three variants were 
nearly close to 1/3 and one variant was elongated along 
the y axis. In the very early stage e < 2 x 10~ 3 , the sys- 
tem deformed elastically with /i e = \l ~ 2. However, in 
the next stage 2 x 10 -3 < e < 0.075, the fraction of the 
favored variant increased up to unity with \x e ~ 0.1 — 0.8. 
This inter-variant transformation occurred without dis- 
location formation. On the return path, the solid was 
composed of the favored variant only with large fi e ~ 7. 
As a a — »• 0, there remained a remnant strain about 0.06. 
However, upon heating to T = 0.1 above the transition, it 
disappeared and the solid again assumed a square shape. 
We note that plastic deformations should occur at high 
strains. In the present simulation, dislocations were in- 
deed proliferated for a a > 0.4 (or e > 0.08) at T = 0.02. 



IV. GLASS FORMATION WITH LARGE 
REPULSIVE IMPURITIES 

In Figs. 7-11, we treat mixtures of elliptic particles and 
large repulsive impurities. With increasing the impu- 
rity composition c, the orientation disorder is enhanced 
and the long wavelength orientation fluctuations are sup- 
pressed, resulting in "orientation-strain glass". Here, 
even for our anisotropic particle systems, we predict the 
nonlinear mechanical behavior studied for strain glass^D. 
The BKT phase disappears with increasing c. 



A. Orientation-strain glass 

Figures. 7-9 are simulation results with a thermostat 
at fixed volume under the periodic boundary condition, 
where T = 0.05 and x = 1-2. The temperature was low- 
ered from a high temperature as in the previous section. 
The size ratio is fixed at o^/ci = 1.2. For c < 0.2, the 
system still forms a single hexagonal crystal with point 
defects at the impurity positions. 
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FIG. 7. Frozen patterns of angles 9i in the range < Qi < n 
(left) and order parameter amplitudes ql (right) with impu- 
rities (black points) at c = 0.05, 0.1 and 0.2, where T = 0.05 
and x = 1-2- 



In Fig. 7, we present snapshots of 9i and qf for three 
compositions as in Fig.l. In the top panel at c = 0.05, 
the impurities induce irregular orientation disorder, but 
not much affect the overall order such that large-scale do- 
mains are still distinct. In the lower panels with c =0.1 
and 0.2, the orientation disorder increases and the do- 
main sizes become finer. For c = 0.2, the system ap- 
proaches orientation glass but with mesoscopically or- 
dered regions still remaining. In Fig. 8, increasing c gives 
rise to suppression of Sq(Ic) at long wavelengths (k < 0.3) 
and (q 2 ) in Eq.(18) at low T. 

Figure 9 displays expanded snapshots of the elliptic 
particles around the impurities. We recognize that the 
alignments are mostly perpendicular to the surface nor- 
mals, analogously to the parallel anchoring of liquid crys- 
tal molecules on the colloid surfaces^. Moreover, we 
notice an apparent tendency of string-like clustering or 
aggregation of the impurities. They tend to be localized 
along the interface regions between different variants, al- 
lowing formation of mesoscopically ordered regions of the 
elliptic particles even for c = 0.2. 

To examine the degree of clustering, we may group the 
impurities into clusters. Let the two impurities i and j 
belong to the same cluster if their distance is shorter than 
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FIG. 9. Expanded snapshots of #,: around large impurities (•) 
for c = 0.1 and 0.2, exhibiting planar anchoring of molecular 
orientations and clustering. 



1.6o"i. Then we obtain the numbers N c \(£) of the I clus- 
ters (those consisting of I impurities), where i — 1, 2, • • • 
and J2e £N c \(l) — N 2 - The probability that one impurity 
belongs to one of the £ clusters is P c \{t) = £N c i(£)/N 2 . 
The average cluster size is defined as 



i cl = j2^dW =J2 £2n ^/ N2 
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FIG. 10. Shape memory effect under uniaxial deformations 
along the y axis with impurities, where T = 0.05, c = 0.2, and 
X = 1.2. Top: e vs a a (left) and (g 2 ) vs a a (right). For 0.12 < 
a a < 0.155, e and (q 2 ) increase steeply. For a a > 0.15, there 
remains only the variant elongated along the y axis. After this 
cycle, the residual strain is 0.04, which vanishes upon heating 
to T = 0.1. Bottom: Snapshots of Oi for a a = 0.12,0.135, 
and 0.15 in the transition region, where large-scale orientation 
fluctuations can be seen but there is no dislocation. 



(33) 



In Fig.7, we have 4i = 1-37, 2.04, and 4.79 for c = 0.05, 
0.1, and 0.2, respectively. 



B. Mechanical properties in glass 

We also observed a shape-memory effect even in orien- 
tation glass, where small disfavored domains were grad- 
ually replaced by favored ones upon stretching. In this 
effect, no dislocation was formed at low T. 

In Fig. 10, at T = 0.05, we increased a a slowly at 
da a /dt — 4 x 10 -6 from up to 0.2, where the variant 



elongated along the y axis becomes increasingly dom- 
inant. We then decreased a a back to at da a /dt = 
—4 x 10~ 6 . Between these two paths, significant differ- 
ences can be seen in the degree of orientation (q 2 ) and 
the strain e. The e is given by Eq.(27) with L y o be- 
ing the initial system length at T = 0.05 and a a = +0. 
On the stretching path, there appear four stress ranges: 
[i e ~ 3 for < a a < 0.05, fi e ~ 0.8 for 0, 05 < a a < 0.12, 
fi e ~ 0.2 for 0.12 < a a < 0.15. Remarkably, the response 
is elastic in the first range and is very large with e in- 
creasing steeply from 0.028 to 0.064 in the third range. 
For a a > 0.15 and on the return path, /i e is of order 



unity and we can see considerable variations in e and 
(g 2 ), where the fractions of the disfavored variants sig- 
nificantly change around the impurities. In contrast, in 
the one-component case in Fig. 6, we have found no such 
changes once a single-variant state is realized. 

In the bottom panels of Fig. 10, we display snapshots of 
0i at four points A, B, C, and D where a a = 0.12, 0.135, 
1.5, and 0, respectively. See the bottom left panel of 
Fig. 7 for the snapshot at the initial time in the same 
run. Between A and B, the orientation and the strain in- 
crease abruptly. In this transition region, we notice emer- 
gence of large-scale orientation fluctuations taking stripe 
shapes and making angles of ±7r/4 with respect to the x 
axis. In stress and thermal cycles in glass, the impurities 
pin the orientation fluctuations in quasi-stationary states 
under very slow time variations of a a and T, yielding the 
history-dependence of the physical quantities. 



C. Positional disorder for <72/ci = 1.4 

So far, the crystal structure has been little affected by 
the orientation fluctuations at <J 2 j<J\ = 1.2 for c = 0.1 
and 0.2. However, if we adopt a larger size ratio and/or 
a larger composition, the positional (structural) disor- 
der is increasingly enhanced, resulting in usual positional 
polycrystal or glass. In our case, the orientation disor- 
der is more enhanced than the positional disorder. This 
is in sharp contrast to liquid crystal systems where the 
nematic order precedes the crystallization. 

In Fig. 11, we set C2/01 = 1.4 and \ = 1-2 to obtain 
polycrystal for c = 0.1 and 0.2. In the left, the orien- 
tation angles 9j are displayed, where there still remains 
noticeable mesoscopic orientation order. In the left, six- 
fold bond-orientation (crystal) angles ctj are displayed, 
where we introduce a ,- for each elliptic particle j in the 
range < ocj < tt/3 by^EH 



2 exp[6i$jk] = Zj exp[6iatj], 



(34) 



fc^bondcd 



Here, Ojk is the angle of the relative position vector 
r jk = r k — r j with respect to the x axis, the particle 
k is within the range |r,-fe| < 1.7a±, and Zj and 6a j are 
the absolute value and the phase angle of the left hand 
side, respectively. For c = 0.1, one large grain is embed- 
ded in a crystal containing many point defects, where 
angle differences are of order 10 — 15 degrees. On the 
other hand, for c = 0.2, many grains appear with much 
larger angle differences. 



ROTATIONAL DYNAMICS 
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FIG. 11. Orientation angle 8i in the range < 9i < -k (left) 
and sixfold bond orientation angle on in Eq.(34) in the range 
< 6i < 7r/3 (right) in polycrystal states for c = 0.1 (top) 
and 0.2 (bottom), where \ — 1-2 and T = 0.05. The size ratio 
is increased to 02/ui = 1.4. 



dependent angle-distribution function defined by 
G (t^) = W E <<%•(* + *<>) -* 3 - («o) -¥>)>, (35) 



Nt 



i<i<jvi 



where the average (• • •) is taken over the initial time to 
and over several runs. Here, G(t,ip) tends to S((p) as 
t — > and is broadened for t > 0. In particular, we treat 
the first two moments G\{t) and G%(t) written as 



2tt 



G\{t) — I dtpG(t,(p)cos(<p), 



/■2-K 

G 2 (t) = I d<pG(t,<p)coB(2(p). 

Jo 



(36) 



(37) 



Since these two functions are unity as t — > 0, we introduce 
two relaxation times, t\ and t 2 , by 



G 1 (r 1 ) = e-\ G 2 {t 2 ) = c- 1 . 



(38) 



These times grow as T is lowered. We plot G\{t) and 
G 2 (t) vs t in Fig.12 and n vs T in Fig.13. 



A. Angle relaxation functions 



B. Turnover motions and configuration changes 



We n ow di scuss the rotation dynamics of the elliptic 
particle d 52 * 53 l In two dimensions, we consider the time- 



As a marked feature, the elliptic particles sometimes 
undergo the turnover motion 8j — > 9j ± it or rij — > —rij 
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FIG. 12. Orientation relaxation functions G\(t) (left) and 
G 2 (t) (right) in Eqs.(36) and (37) for c = and X = 0.6 (top) 
and for c = 0.2, cra/ci = 1.2, and x — 1-2 (bottom). Relax- 
ations are slowed down as T is lowered. The Gi(t) decays due 
to turnover motions, while G2(t) due to configuration changes. 
For c > 0, G2(t) tends to a finite constant as i — > oo. 
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FIG. 13. Orientation relaxation time n from Gi(t) for (a) 
c = and \ — 0-6, (b) c = 0.1 and \ = 0.6, (c) c = and 
X = 1.2, and (d) c = 0.1 and \ = 1-2- It is the time scale of 
successive turnover motions. For (a) and (c), n grows steeply 
in the BKT phase (T 2 < T < Ti). For strain glass (b) and (d), 
the BKT phase is nonexistent and ti grows as T is lowered. 



taking place in a microscopic time (~ lp^l. In terms of 
the orientation vector rij , we also have 

G?i(t)= ^ (n^t + to)-^-^))/^, (39) 

1<J<»1 

so that the times between successive turnovers of an el- 
liptic particle are of order t\ in Eq.(38). On the average 



over all the elliptic particles, the turnover motions give 
rise to a peak in G(t, ip) at ip = it. 

In our simulation, it is nearly of the Gaussian form for 
t <C T\ in the range \ip — tt| < 1 as 



G(t,^) = 



A(t) 

2lT(T 



exp 



&-*? 



2ct 2 



(40) 



where the variance a is a constant about 0.45 in the 
present case. The integral of this Gaussian peak is equal 
to the coefficient A(t), so A(t) has the meaning of the 
turnover probability per elliptic particle in the time in- 
terval [0, £]. In terms of ti, we find the linear growth, 



A(t) * Cit/n, 



(41) 



in the early time range t <C T\. In our system C\ = 0.5. 
On the other hand, G2{t) is unchanged by the instanta- 
neous turnover motions, so it relaxes due to the orien- 
tational configuration changes involving the surrounding 
particles. We found the inequality r 2 > t\ at any T and 
c in our simulation. 

In Fig. 12, both Gi(t) and G 2 (£) relax considerably in 
the early time region t < 2 due to the thermal rapid mo- 
tions of the orientations without configuration changes. 
For t > 2 the fitting G x {t) ~ exp[-(t/ri) /3 ] fairly holds, 
where (3 decreases from unity to about 0.5 as T is low- 
ered. Furthermore, for c > 0, G^{t) tends to a nonva- 
nishing positive constant fa at large p2. In our case, 
this plateau appears because the anchoring of the elliptic 
particles around the impurities becomes nearly perma- 
nent at low T. Thus we found that the plateau height / 2 
increases with lowering T and with increasing c. 

In Fig. 13, the two curves for c = indicate that t\ is 
short (< 10) for T > T\, increases steeply in the BKT 
region T% < T < Ti, and grows further in the ordered 
region T < T 2 in the thermal activation form, 



ti ~ exp(r /T) (T<T 2 ). 



(42) 



We have T ~ 0.40 at \ = °- 6 on curve (a) and T ~ 1.2 
at x — 1-2 on curve (c). In addition, n ~ r 2 for T > Ti 
but t 2 /ti > 1 for T < T 2 . In fact, for x = 0.6, t 2 /ti is 
about 10 2 at T = 0.07 and is about 10 3 at T = 0.06. 

For c > 0, the turnover motions still occur with t\ < 
t-2- However, in Fig. 13, the relaxation behavior for c > 
is very different from that for c = 0. In the disordered 
phase with T > T\ , t\ for c > is longer than n for c = 
due to the impurity pinning. For T < T 2 , on the contrary, 
n for c > is shorter than n for c = 0. That is, the 
turnover motions are more frequent in orientation glass 
with c > than in the oricntationally ordered phase with 
c = 0, as ought to be the case. Above T 2 , the impurity 
anchoring gradually becomes transient. 

It is worth noting that Chong et aP^ studied the 
orientation dynamics of a glass-forming binary mix- 
ture of dumbbells using the angle relaxation functions 
C e {t) = J2j( p e{nj{t + t) -rijita^/N in three dimen- 
sional molecular dynamics simulation, where Pi is the 
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FIG. 15. Expanded snapshots of elliptic particles and small 
attractive impurities (black points) in the boxes in Fig. 14, 
where homeotropic anchoring and impurity clustering are 
marked. Colors of the elliptic particles represent the ori- 
entation angles 8j (left) and the sixfold crystal angles otj 
(right) from the same data. The domains (left) are finer than 
the grains (right). The cooling rate from T = 1 to 0.1 is 
dT/dt = -1.8 x 10~ 5 . 



FIG. 14. Snapshots of orientation angles 0i (left) and sixfold 
crystal angles o^ with addition of small attractive impurities 
with c = 0.1 (top), 0.2 (middle), and 0.3 (bottom), where 
T = 0.1, x = 1-2, 0-2/0-1 = 0.6, and ( = 2. The orientation 
disorder is stronger than the positional disorder. The cooling 
rate from T = 1 to 0.1 is dT/dt = -1.8 x 10" 5 . 



Legendre polynomial of order £ and rij is the orientation 
vector of particle j. The relaxations of C\{f) and C-^ii) 
for small dumbbell anisotropy in their paper closely re- 
semble those of G\(t) and Gi{t) for c = 0.2 in Fig. 12. 



VI. GLASS FORMATION WITH SMALL 
ATTRACTIVE IMPURITIES 

In this section, we further treat another intriguing case 
of small attractive impurities with 02/ a\ = 0.6 in Eq.(3) 
and with £ = 2 in Eq.(6). Such small impurities tend to 
be expelled from the ordered domains of the host parti- 
cles. We shall see that they form clusters. 



A. Orientational disorder and positional disorder 

Though not shown in this paper, we performed simula- 
tion runs for small repulsive impurities with o~i/ <j\ = 0.6 
and C, = 0, where most of the impurity aggregates are 
stringlike and the anchoring of the elliptic particles is 
planar. However, if the anisotropy strength of attraction 
C, is increased at fixed \, the aggregates becomes increas- 
ingly compact. For £ > 1, the aggregates can "solvate" 
several elliptic particles in the homeotropic alignment^. 
With further increasing £, even a single impurity creates 
a solvation shell composed of several elliptic particles like 
a small metallic ion in water. 

In Fig. 14, we show snapshots of 9j and otj of all the 
particles, where T = 0.1, x = 1-2, 02/ci = 0.6, and 
C = 2. Here, we set dT/dt = -1.8 x 10~ 5 . For c = 0.1, 
the system is still in a single crystal state, but the orien- 
tational domain structure induces large-scale elastic de- 
formations, leading to close resemblance of the patterns 
of 9j and oij. For c = 0.2, the orientational domains are 
much finer and a polycrystal state is realized with larger 
grains (> 10). For c = 0.3, the orientation order is much 
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FIG. 16. Snapshots of orientation angles 9i (left) and sixfold 
crystal angles Ui (right) with addition of small attractive im- 
purities with c = 0.1. The parameter values are common to 
those in the top panel of Fig. 14, but the cooling rate from 
T = 1 to 0.1 is dT/dt = -9 x 10" 3 . Here, the degree of 
clustering is weaker, resulting in a polycrystal state. 
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FIG. 17. Probability P cl (£) = £N cl (£)/N 2 of an impurity 
belonging to £ clusters, where N c i(£) is the cluster number 
composed of £ impurities. Here small attractive impurities 
are considered for slow cooling in Fig. 14 and fast cooling in 
Fig. 16. Distribution is broader for slower cooling. 



more suppressed and a positional glass state is realized 
with mesoscopic heterogeneities still remaining. 

In Fig. 15, we display expanded snapshots of 8j (left) 
and otj in Eq.(34) (right) in the box regions in Fig. 14. 
The alignments of the elliptic particles around the impu- 
rities are mostly parallel to the surface normals. This is 
analogous to the homeotropic anchoring of liquid crystal 
molecules on the colloid surfaces^. We notice a tendency 
of clustering or aggregation of the impurities. Compar- 
ing the left and right panels, we recognize that the inter- 
faces are finer than the grain boundaries. That is, the 
interfaces can be seen both on the grain boundaries and 
within the grains. The impurities tend to be localized on 
the interface regions between different variants. 



B. Cooling-rate dependent clustering of impurities 

The degree of impurity clustering should be decreased 
with increasing the cooling rate dT/dt for long diffusion 
times of impurities. In Fig. 16, dT/dt is —9 x 10~ 3 and is 
500 times faster than in Fig. 14, where the other param- 
eters are common. We give snapshots of 9j and ctj at 
c = 0.1, where the clustering can be more evidently seen 
than for c = 0.2 and 0.3. While a single crystal has been 
realized in the top panel of Fig. 14, a polycrystal state is 
realized with large angle differences in Fig. 16. 

Let the two small impurities i and j belong to the same 
cluster if their distance is shorter than \.1o\. Then we 
obtain the number N c \(£) of clusters composed of £ im- 
purities. In Fig. 17, we show the cluster size distribution 
P cl (£) = £N cl (£)/N 2 {I = 1,2,- ••) for the_ examples in 
Figs. 14 and 16. The average cluster size £ c \ in Eq.(33) 
increases with c as 2.45 for c = 0.1, 3.43 for c = 0.2, and 



4.61 for c = 0.3 under the slow cooling in Fig. 14, while 
t c \ = 1.62 for c = 0.1 under the fast cooling in Fig. 16. 

It is knowrpS that water becomes glass at low T 
with addition of a considerable amount of LiCl, where 
small hydrophilic Li + and Cl _ ions solvate several water 
molecules via the strong ion-dipole interaction. The re- 
sultant orientation anchoring of water molecules should 
even prevent formation of the crystal order at high salt 
concentrations, resulting in the observed positional glass. 
It is natural that the cooling rate influences the degrees 
of ion clustering and vitrification. 



VII. SUMMARY AND REMARKS 

We have presented an angle-dependent Lennard- Jones 
potential for elliptic particles and impurities, which 
depends on the orientation angles of the interacting 
particles. Using this potential, we have performed 
simulation of 4096 particles on very long time scales 
(~ 10 5 To) in two dimensions. Our main results are as 
follows. 

(i) In Sec. II, we have presented our model potential, 
where the anisotropy strengths are characterized by x 
for the repulsive part in Eq.(5) and C, for the attractive 
part in Eq.(6). The aspect ratio of the elliptic particles 
is given by ag/a s = (1 + 2x) 1 ^ 6 - In this paper, \ is °f 
order unity, so we have assumed weak particle anisotropy 
to find crystallization at a high temperature above the 
orientation transition. 

(ii)In Sec. Ill, we have presented simulation results for 
one-component systems of elliptic particles by changing 
the temperature T to produce Figs. 1-6. The domain pat- 
terns in Fig.l at low T are those observed on hexagonal 
planes. In our case, the Berezinskii-Kosterlitz-Thouless 
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phase32ESI [ s realized in a temperature window, where 
the orientation fluctuations are much enhanced at long 
wavelengths as indicated by the structure factor Sq(Iz) 
in Fig. 3. We have shown thermal hysteresis in Fig. 4, 
singular behaviors of the shear modulus and the specific 
heat in Fig. 5, and a shape-memory effect in Fig. 6. 
(iii) In Sec. IV, we have examined orientation-strain 
glass of elliptic particles and large repulsive impurities 
with the size ratio o^/ci = 1.2 in Figs. 7-10. The 
orientations of the elliptic particles are pinned at the 
impurity surfaces in the planar alignment in Fig. 9. The 
shape- memory effect in strain glass is marked in Fig. 10. 
Positional disorder also emerges for o<zl<j\ = 1.4 in 
Fig. 11. 

(iv) In Sec.V, we have studied the rotational dynamics 
of the elliptic particles. In Fig. 12, G\(t) decays due to 
the turnover motions of the elliptic particles, while G2(t) 
decays due to the configuration changes. In Fig. 13, the 
turnover relaxation time t\ grows at low T and behaves 
differently with and without impurities, 
(v) In Sec. VI, we have examined the effect of small 
attractive impurities on the orientation disorder and 
the positional disorder in Fig. 14. The impurity effect 
is stronger on the former than on the the latter. The 
elliptic particles are homeotropically anchored at the 
impurity surfaces in Fig. 15. The clustering of impurities 
is suppressed for rapid cooling as in Figs. 16 and 17. 



We further make critical remarks as follows: 

(1) In our simulation, we used a Nose-Hoover thermostat 
(NHT f® in all the figures and a NHT and a Parrinello- 
Rahman barostatpEI in Figs. 4-6, and 10. In future work, 
we should examine the coupled dynamics of the trans- 
lational and orientational degrees of freedorrP without 
thermostats and barostats in the system interior. 

(2) There has been no systematic measurement of the 
mechanical properties of orientationally ordered, multi- 
variant crystal and orientation glass. Such experimental 
results could be compared with those from shape-memory 
alloys^-U Weak elasticity was observed in orientation- 
ally disordered solids above the transition (called "plastic 
solids" ) in creep experiments^. Also, as far as the authors 
are aware, there has been no experimental information of 
the impurity clustering in any physical systems exhibit- 
ing mesoscopic heterogeneities. 

(3) In this paper, the particle anisotropy is not large, 
which favors formation of crystal order. For large 
anisotropy, liquid crystal phases should appea r 41 * 42 ^, 



where the impurity effect is of great interest. As sug- 
gested by the experiment^, addition of a considerable 
amount of impurities leads to the orientation order only 
on mesoscopic scales in liquid crystal phases. In such 
states, we expect large response to applied electric field. 

(4) In Figs.l and 4, our system undergoes a structural 
phase transition gradually in a narrow temperature win- 
dow even for the one-component case. In our model, a 
gradual phase transition still occurs in the stress-free con- 
dition without impurities. However, we also stopped the 
cooling and waited for a long time (3> 10 4 ) at T — 0.06 on 
the stress-free cooling path in Fig. 4; then, we observed a 
transition to the ordered single- variant phase (not shown 
in this paper). Thus, in future work, we need to calculate 
the Gibbs or Helmholtz free energy to decide whether the 
system is in equilibrium or in a metastable state. 

(5) As well as the orientation fluctuations, the displace- 
ment fluctuations are also enhanced around the orienta- 
tion transitiomas indicated by Fig. 5 and by the previous 
experiment d 1 | 4 | 5 l In addition, according to Cowley's clas- 
sification of elastic instabilities^, our phase transitions 
belong to type-I instabilities where acoustic modes be- 
come soft in particular wave vector directions. 

(6) The disordering effect induced by impurities prevents 
a sharp transition^. Thus there is no sharp phase bound- 
ary between the high-temperature orientationally disor- 
dered phase and the low-temperature orientation-strain 
glass phase. These two phases change over gradually with 
varying T as in the cases of positional glass transitions. 

(7) Ding et alP^ numerically studied the superelastic- 
ity, whic h arise s from a stress-induced martensitic phase 
transitio n 1 10 * 11 !. We also realized this phenomenon for 
anisotropic particles, which will be reported elsewhere. 

(8) We will also report three-dimensional simulation on 
mixtures of spheroidal particles and spherical ones with- 
out and with the dipolar interaction. We shall see finely 
divided domains produced by impurities and large re- 
sponses to applied strain and electric field. 
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